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Abstract — We extend the Wiener-Khinchin theorem to non- 
wide sense stationary (WSS) random processes, i.e. we prove that, 
under certain assumptions, the power spectral density (PSD) of 
any random process is equal to the Fourier transform of the 
time-averaged autocorrelation function. We use the theorem to 
show that bandlimitedness of the PSD implies bandlimitedness of 
the generalized-PSD for a certain class of non-WSS signals. This 
fact allows us to apply the Nyquist criterion derived by Gardner 
for the generalized-PSD. 

Index Terms — Non-wide Sense Stationary Processes, Power 
Spectral Density, Subsampling, Wiener-Khinchin Theorem, Ban- 
dlimlted 



I. Introduction 

The Power Spectral Density (PSD) defined in dU of a 
random process is the expected value of its normalized pe- 
riodogram, with the duration over which the periodogram 
is computed approaching infinity. For wide sense stationary 
(WSS) processes, the Wiener-Khinchin theorem [6] shows that 
the PSD is equal to the Fourier transform of the autocorrelation 
function(treated as a function of the delay). Wiener-Khinchin 
theorem is a very fundamental result because it can be used 
to perform spectral analysis of WSS random processes whose 
Fourier transform may not exist. 

In this work we give a complete rigorous proof of the 
nonstationary analog of the Wiener-Khinchin theorem, i.e. we 
show that under certain assumptions, the PSD of a random 
process, defined in (1), is equal to the Fourier transform of 
the time-averaged autocorrelation function. 

While similar ideas have been introduced in some textbooks 
on random processes[l][4][10][ll], (with the aim of gener- 
alizing the stationary case), with the exception of [11], the 
results are mostly incomplete. We discuss this in section 2.1 
and give details in [12]. Also, the result of [1 1] uses a different 
sufficient condition than ours. Note that the work of [2] has a 
very related title, but a completely different contribution from 
the current work. It proves Wiener-Khinchin theorem for a 
generalization of autocorrelation. 

There are other commonly used definitions of the PSD 
for non-WSS processes, e.g. the generalized-PSD [3], K(u,v), 
which is the 2D Fourier transform of the autocorrelation or the 
evolutionary spectral density function [5]. In [8], K{u,\') has 
been shown to be equal to the covariance between the Fourier 
transform coefficients of the signal at u and at v. 

Our result is important because one can now do 1- 
dimensional spectral analysis of certain types of non-WSS 
signals using the PSD. The generalized-PSD, K(u, v), defines 



a 2D Fourier transform which is much more expensive to 
estimate (as explained in section 3) and more difficult to 
interpret. We prove in section 3 that at least for certain 
types of non-WSS signals, the bandwidth computed using 
the PSD can be used for subsampling the signal. This is 
done by using our result to show that bandlimitedness of the 
PSD implies bandlimitedness of the generalized-PSD for a 
certain class of non-WSS signals and then using Gardner's 
result [3] for subsampling. One motivating application for the 
above is analyzing piecewise stationary signals for which the 
boundaries between pieces are not known. It is computation- 
ally more efficient if one can first uniformly subsample the 
signal to perform dimension reduction (by using the bandwidth 
computed from the estimated PSD) and then use existing 
techniques to find the piece boundaries or to perform inference 
tasks such as signal classification. 

A common practical approach for computing the PSD of 
a WSS process is to break up a single long sequence into 
pieces and to use each piece as a different realization of the 
process. But this cannot be done for a non-WSS process. 
Multiple realizations are required before either the autocor- 
relation or the PSD can be estimated and these are often 
difficult to obtain in practice. One application where multiple 
realizations are available is to analyze time sequences of 
spatially non-WSS signals, which are temporally independent 
and identically distributed (i.i.d) or temporally stationary and 
ergodic. For example, if the sequence is a time sequence of 
contour deformations or a time sequence of images (2D spatial 
signals) when temporal i.i.d-ness or stationarity is a valid 
assumption. The contour deformation "signal" at a given time 
is a ID function of contour arclength and is often spatially 
piecewise stationary, e.g. very often one region of the contour 
deforms much more than the others. Temporal stationarity 
is a valid assumption in many practical applications such as 
when analyzing human body contours for gait recognition or 
analyzing brain tumor contour deformations and thus the time 
sequence can be used to compute the expectations. 

Paper Organization: In Sec. 2, we give the complete 
proof of Wiener-Khinchin for non-WSS processes. In Sec. 3, 
sampling theory for a particular class of non-WSS processes 
is discussed. Application to time sequences of spatially non- 
WSS signals is shown in Sec. 4. Conclusions are given in Sec. 
5. 
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II. NoN-WSS WiENER-KHiNCfflN THEOREM Sincc R^iT) is integrable, we can use dominated convergence 

J, , / N • I 1 theorem [7, Chap 4] to move the limit inside the outer integral 

The PSD or any random process, x(f), is denned as • , ^ r r,,, , • , 

in the first term or (qJ. Thus, this term becomes 



(8) 



E[\Xt{oj)\^] 1 
^■'■(^^ = ^™o 2f ^ Xt(oj)^ x(t)e-J dt (1) ]im[±, RMJ2-T)dt2gT{T)]e-^'''d 

' J-oo r^"" J-T 

E[.] denotes expectation w.r.t. the pdf of the process x{t). Let C'^ ^. 1 f f )dt ] ^^'^^d 

■/^xCfi, fi) - £'[^(ri)^(f2)] be the autocorrelation of x(f). Assume ~ J-co 2r ^ ^ 

r r £[|x(fi)x(f2)|](ifiiif2 < "3 (2) 

J-T J-T 

for any finite 70. Because of (2), Fubini's theorem [7, Chap 

12] can be applied to move the expectation inside the integral ,. 1 T^' ' ,n ^ \ -wjti i i /r.^ 

. . . . . , • . • X. < lim — \RAt2,t2-T)e ' \dt2dT (9) 

and to change the order of the integrals in (1). Also, assume r^oo 2T Jj^q Jt=-T 

that the "maximum absolute autocorrelation function", p^.^^^ q, we know \RM, ti - t)| < R":{r). Hence 

R":(T)^mp\RAt,t-T)\ (3) 



Now, for the second term in (|5]), we know 

Um|-^ r f RAtiJ2-T)e-i"'dt2dT\ 
■-IT rT-T 



lim|-^r r R,{t2,t2-T)e-i''-dt2dr\<\im { t^^Wt 
r^oo 2T J^^o Jt2=-T Jt=o 2r 

is integrable. Moving the expectation inside the integrals in (10) 
O and defining r ^ - ?i , we get ^ ^^^^^ if < r < 2r 

^ ^ [0 otherwise 

^-^■("^ = JT. 4^1).-" ^'^i £ ^j^^^ j.^^^^ j^2r ^^^^^^^^ ^ j.^^^^ 0,(r)/^?(r)^.. 

1 We can apply dominated convergence theorem again since 

= lim — I I R,(t2,t2-T)e-'"^dT dt2 (4) \cj)j(T)R";(T)\ < R"'(t) which is integrable. Thus, this implies 

Changing the integral order and the corresponding limits of lim I 4>T(T)R'"(T)dT - I lim [4>t (t)R'^ (T)]dT - (11) 

7"— >oo I ' I 7" — >oo 

integration, in a fashion similar to the proof of the stationary 

case [6], Therefore, ( fTOb forces 

5.(c.) = lim -L f r' ;^,(,,, - r)e-'-^dt2dr + « < lim J ' RAt2, t2 - r)e-^-^ dt2dr\ < 

7-^0= 2T Jr^o J,,=r-T 2T T-T 

lim -L f r ' RAt2, t2 - r)e-"^^dt2dr ^ 2^ L f r " -)^-^'"^*2^r = (12) 

Y n2T pT Thus, the second term of (|5]l is 0. In an analogous fashion, we 

~ T^ 27 J 27- J T ^ ~ ^''^ ^'^''dt2dT - can show that the third term is also 0. Thus, we finally get 

lim — r r R,{t2,t2-T)e-j"'dt2dT- 5,-(w) = R,{T)e--""'dT, where (13) 

lim -L f f RAt2, t2 - T)e-'-^dt2d. (5) ^-^-W - ^ ^^(^2, ^2 - t)^^^^ (14) 

This is the Wiener-Khinchin result for any general random 

Define (t) ~ 1^ < t < 2r process, i.e. 

1^0 otherwise Theorem 1 {Wiener-Khinchin Theorem for Non-WSS Processes): 

^, , „ r , If ^ holds, and if 7?™(t) defined in OJ is integrable, then 

Then the first term of (qJi becomes u u ■ .u ncr> j ^ j • m • T . r- 

( 1131 ) holds. I.e. the PSD, defined in (UJ, is equal to Fourier 

1 . transform of the "averaged autocorrelation function", defined 

^ 1-2T i-T " "^^""^'^^^^^ in m. 

^ fT(T)dT (6) ^ Discussion of related results 

Related ideas are introduced in several textbooks[l] [4] 
where /^(t) = ^ J^j. R.x{t2,t2 - T)dt2 grir) e-j'^\ It is easy to [10][11] , but except for [11], the results in the rest of them 
see are either different from ours or incomplete. 

2J. [1] does not justify why the limit can be moved inside the 

l/r(T)| < ^R'^(T)gT(T) = R'^(T)gT(T) < R'^(t) (!) integral in equation 7.38. [4] has a similar problem and also 

his result says that S{f) is equal to the time average of 
'A sufficient condition for (2) is: ffiere exists a C < c» such that rz) < *^ instantaneous PSD. The condition of [10] given on page 
C almost everywhere (except on a set of measure zero) 179 is not sufficient either. Fourier transform of the averaged 
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autocorrelation function may not exist even if his condition 
is satisfied. Peeble's result[ll] gives a different sufficient 
condition from our result (He assumes absolute integrability 
of PSD). We discuss the results of the above books in detail 
in [12]. 

III. SUBSAMPLING A ClASS OF NoN-WSS PROCESSES 

We use Gardner's result [3] for subsampling nonstationary 
random signals along with the PSD, instead of the generalized- 
PSD, since the PSD is much less expensive to compute. For a 
T-length signal for which P realizations are available, the PSD 
defined in (1) can be estimated in 0{PT log T) time (need to 
estimate P T-length Fourier transform and average their square 
magnitudes). Generalized-PSD computes RxitiJi) using P 
realizations (takes 0{PT^) time), followed by computing its 
2D Fourier transform (takes 0{T^ log T) time), i.e. it requires 
0(PT^ + log T) time. Clearly O(PTlogT) < O(PT-). 

Gardner's result says that if the generalized PSD of x, 
K^iu, v), is bandlimited in both dimensions, x can be recon- 
structed exactly (in the mean square sense), from its uniformly- 
spaced samples taken at a rate that is higher than twice the 
maximum bandwidth in either dimension. To use this result 
with the PSD, we need to show that bandlimitedness of the 
PSD implies bandlimitedness of the generalized-PSD. The 
most general case for which this can be done will be studied in 
future. We show it for the following class of non-WSS random 
signals, which can be used to model many commonly occur- 
ring random processes including many piecewise stationary 
ones. This is one of the four classes of nonstationary processes 
described in [9]. 

Definition 1 ( Class NSl ): A random signal, x{t), belongs 
to the class NSl if it can be represented as the output of 
"nonstationary white Gaussian noise", w(f), passed through 
a stable linear time invariant (LTI) system, denoted h{t), i.e. 
x{t) - w{t) * h(t), where ★ denotes convolution; w(f) is a 
Gaussian process with E[w{t)] - 0, Vf, Rn{t, t') - (T^^,{t)6{t-t') 
and cr~^,,„ax - max, crj^,(0 < and hit) satisfies j^^ \h{t)\dt < 

oo. 

We show below that for this class of signals, Rx{t) - ct^Mt) ★ 
h{—T) and that Theorem[T]can be applied to show that 5^(w) = 
0-2 |//((j>)p where H(lo) denotes the Fourier transform of h{t). 
Using the definition of NS 1 signals, Rx{t) is written as 

f^T r^oo /^oo 
T^oo 21 J_J 

w(ti - T - T2y\dT2dr\dt\ 

= lim — h{Tl)h{Tl - T)(rl,{tl - Tl)dTldtl (15) 

By changing the order of integration and defining t - t\-T\, 
we get 



RAt) = lim 



fhf" 

J-oo ^-f J-r-T, 



Crlmt h{T,)h{T,-T)dTi{\6) 



The integrand, i/'(t,Ti) = ^ J_^^^^crl{t)dt h{Ti)h{Ti - t) < 
cri,^maxh{T\ - T)h(r\). Since h{T\ - t)/z(ti) is absolutely inte- 
grable w.r.t. ti (follows from Cauchy-Schwartz inequality and 
the fact that a stable h(t) implies that h^(t)dt < oo), we 



can apply the dominated convergence theorem [7] to move 
the limit inside to get 

lim{— Crl(t)dt} h{Tl)h(Tl-T)dTl(n) 

^r^oo 21 J-T-Ti 
By using an argument similar to that used in (fT2l i and the 
paragraph below it, we can replace lim^^oo S-t-t ^'^ 
Hmj-^oo Tr J_7-[-]- Thus, we get 

(^lMTl)hiTl - T)dTl = alh(T) * /i(-t),(18) 

CO 

1 

where al = Um — crl(t)dt (19) 
r-ioo 21 J_j 

By taking the Fourier transform of both sides of ( fTSl l. and 
using Theorem [T] for the left hand side, we get that 



(20) 



Theorem [T] can be applied because (a) R'^ij) - sup, \Rx{t, t - 
t)| < cr^, „jgJh{T) -k h{-T)\ which is absolutely integrable (since 
h(T) ★ h(-T) is stable) and (b) the sufficient condition given 
in footnote 1 for ^ to hold is satisfied since R\x\{t\,t2) < 
sup, Ri4t,t) = sup,7;^(f,f) < 0-2 ,„^, £^/j2(ti)c/ti ^ C. 

From ( l20b . for signals belonging to NSl, if 5i(w) is 
bandlimited, it implies that H{cl)) is bandlimited. This in 
turn impHes that Kx(ii,v) - J J Rx(ti,t2)e'^''"''''''-''dtidt2 - 
H(u)H*{v)Kk{u, v) is also bandlimited with the same band- 
width. Since Kx{u, v) is bandlimited, Gardner's result [3] 
applies. Thus we have the following corollary. 

Corollary 1: For random signals x{t) belonging to the class 
NSl, if S x{(^) is bandUmited with bandwidth B, then x{t) 
admits the following mean-square equivalent "sample repre- 
sentation" 



n{t-nT,) 



E[(x{t) - 2 x{nT,) ^^J^^^ f] = 0, V f 



(21) 



if Ts < l/(2B), i.e. x{t) can be reconstructed "exactly" (in the 
mean square sense) from its samples taken less than 1/{2B) 
interval apart. 

IV. Application to Time Sequences of Spatially Non-WSS 
Signals 

To compute the PSD and its bandwidth for non-WSS 
processes, multiple realizations of the process are required. 
One application obtaining these is to deal with a time sequence 
of spatial signals. We can compute an estimate of the spatial 
PSD using the signals from a temporally i.i.d or a temporally 
stationary and ergodic sequence as the multiple realizations. 
In this section, we demonstrate Theorem 1 and Corollary 1. 

In all simulations, we use discrete time signals. The spatial 
index is denoted by k and the time index by n, i.e. x„[k], as 
a function of k, is a spatial signal. Averaging over time, n, is 
used to compute expectations. For FiglT] we generate x„[k] by 
passing temporally i.i.d. spatially nonstationary white noise, 
w„[k], with o-^,.„[l : A^] = [1^/3,0.12/^/3] (where cj denotes 
a vector of c's of length j), through a stable infinite impulse 
response (IIR) filter, h[k], with transfer function H{oj) = 1/(1 - 
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oie^j" -a2e^^-'") and ai - 0.8, 02 = 0.1. This is done for each 
n, i.e. x„[k] = w„[A:] ★ V«. The signals, x„[A:], belong to 
the class NS 1 and so they satisfy the assumptions of Theorem 
1 (as explained in Section 3). We show the verification of 
Theorem 1 in Fig. [T] The RHS of (fT3T l (Fourier transform 
of averaged autocorrelation function) is plotted as a blue '- 
o' line. Autocorrelation and its spatial average are computed 
using an = 500 length signal (to approximate 00). We also 
plot the A^-length PSD, i.e. the expectation (here average over 
time) of the normalized periodogram of an A^-length signal, 
for increasing values. 

In Fig.|2] we demonstrate Corollary 1 for the same data, but 
now using the knowledge of cr^_„[A:] and h[k]. The theoretically 
computed averaged autocorrelation, crl,h{T) ★ h{—T), and its 
Fourier transform, cr~^\H{w)f', are plotted as a blue '-o' line 
in the top two figures. The numerically computed averaged 
autocorrelation and A^-PSD are also plotted for increasing in 
the same two figures. Next, we use the A^-PSD, for = 500, to 
compute the 90%-bandwidth (point beyond which the residual 
PSD sum is less than 10% of the total PSD sum) and decimate 
all x„'s at twice this rate. The computed MSE between x„ and 
its reconstruction, i„, is 8.05%. This MSE approaches 10% as 
A^ — » 00. Similar results are obtained when we simulated x„[k] 
from a temporally stationary process instead of a temporally 
i.i.d. one. This is done by generating each w„[A:] from a 
temporal AR-1 process and passing it through the same IIR h 
as in Fig. |2l 

V. Conclusions 

We have proved the Wiener-Khinchin result for non-WSS 
processes. This has been combined with Gardner's result [3] 
to prove that Nyquist's criterion can be used to subsample 
a certain class of PSD-bandlimited nonstationary signals. 
Application of these two results to subsampling a simulated 
time sequence of spatially non-WSS signals is shown. Future 
work includes proving a general Nyquist-type result for all 
non-WSS signals that satisfy the assumptions of Theorem [1] 



N-PSD and FT of averaged autocorrelation 



Averaged autocorrelation 




Fig. 1. Demonstrating Theorem 1 for a nonstationary white noise 
passed through an IIR filter (explained in the text). The blue '-o' line 
is the Fourier transform of averaged autocorrelation (RHS of ( 113b ) 
with the averaged autocorrelation approximated using A' = 500-length 
signals. The other five lines are the A',-PSD's with A', increasing from 
A'l = 400 to A'5 = 500. As can be seen, this approaches RHS as A' 
increases (here since RHS is approximated using 500 points, LHS is 
exactly equal to RHS for A'5 = 500). 




50 100 150 200 250 

Fig. 2. Demonstrating parts of Corollary 1. These figures are for 
temporally i.i.d spatially nonstationary white noise passing through 
a stable IIR AR-2 filter with H(oj) = 1/(1 - ciie-i'^ - ciie-J^'"). The 
first figure verifies l llSt . The middle figure demonstrates l l20t and in 
fact also Theorem 1. In the bottom one, one spatial signal and its 
reconstruction by subsampling using a 90%-bandwidth (computed 
using the N=500-PSD) and reconstructing the signal is shown.The 
mean of the square of the residual error over all sequences is 8.05% 
(will approach 10% if the A' in the PSD is increased to 00). 
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